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Abstract — We propose to apply the Back and Forth Nudging 
(BFN) method used for geophysical data assimilations [1] to 
estimate the initial state of a quantum system. We consider a 
cloud of atoms interacting with a magnetic field while a single 
observable is being continuously measured over time using 
homodyne detection. The BFN method relies on designing an 
observer forward and backwards in time. The state of the BFN 
observer is continuously updated by the measured data and 
tends to converge to the systems state. The proposed estimator 
seems to be globally asymptotically convergent when the system 
is observable. A detailed convergence proof and simulations are 
given in the 2-level case. A discussion on the extension of the 
algorithm to the multilevel case is also presented. 

I. Introduction 

Estimating the state of a quantum system is a fundamental 
problem of great interest in quantum control. Amongst a 
variety of applications, it is essential to verify the efficiency 
of a quantum state preparation protocol [2], [3]. For this 
reason, it is interesting to avoid the usual quantum state 
tomography scheme which involves doing the experiment 
many times and performing a strong projective measurement 
of a new observable at each preparation. Indeed, since many 
realizations of the preparation protocol are necessary to 
obtain one state estimation, the fidelity of the preparation 
protocol is averaged out over all these realizations. A new 
approach overcoming this problem was proposed and verified 
experimentally in [4] where a controlled evolution is applied 
to an ensemble average while an observable is continuously 
measured. A Baysesian filter is then used to reconstruct the 
quantum state from the measurement record. In this paper, 
we consider a similar setting to the one in [4], [5], [6]. We 
propose a new approach inspired of the BFN method used in 
geophysical data assimilation [1] to reconstruct the state of 
the system from the measured data. We design an observer, 
which is an estimation of the quantum systems state and feed 
in the data continuously until the observer converges to the 
systems state. A similar proposal was outlined in [7]. How- 
ever, the method we propose has the advantage of extending 
naturally to a multidimensional case and makes more use 
of the specific dynamics which the system undergoes. We 
guess this can strongly reduce the computation time of the 
estimation and increase it's robustness. 
In section O we detail the problem settings and give the 
dynamics equations of the BFN observer Simulation plots 
are then presented in section |ll| to demonstrate the efficiency 
of the state reconstruction protocol. A detailed convergence 



proof of the observer is given in section |iy[ Finally, in 
section ^ we discuss the extension of this algorithm to the 
multilevel case. 

II. The problem setting 

We consider the experimental setting introduced in [4] . To 
simplify the theoretical study, we suppose that the system 
is a spin ^ system (instead of a system of total angular 
momentum equal to 3 or 4 as considered in [4]). It interacts 
with a magnetic field in the x-y plane: the control, and a 
probe. Homodyne detection of the probe, as explained in [5] 
enables a weak continuous measurement of the spin system. 
We suppose that all the parameters involved in the dynamics 
are known. The dynamics of the spin ensemble is described 
by the master equation: 



pit) = -i[B4t)a^ + By{t)ay,p{t)] 

+ r{a,p{t)a, - pit)) 
y{t) = Tr(<7.p(f)) 



(1) 



[.,.] is the Lie Bracket operator and Tr(.) is the trace 
operator We take h ^ 1. p{t) is the density matrix of the 
average ensemble at time t. It is a 2 x 2 positive Hermitian 
matrix of trace 1 . y is the measurement, a^ , <Jy and a^ are 
the standard Pauli matrices. F > gives the strength of the 
coupling between the probe and the system. The aim is, from 
a set of data {y{t)\t e [0, T]}, to estimate the initial state of 
the system p(0) which can be any pure or mixed state. 
In order to do so we use Luenberger observers based on the 
back and forth nudging (BFN) method [1], [7]. The designed 
observer was first introduced, without BFN and for a non 
dissipative system, in [8]. 

The idea is to design an observer on system (jTl) and another 
on the same system but by changing t ^ T — t. And doing 
this iteration n times. This is equivalent to supposing that 
the system has a periodic dynamics of period 2T which is 
symmetric with respect to t = T and that we are measuring 
the system over a time interval [0, 2?iT]. This gives more 
time for the observer to converge with a small gain and with 
minimal amount of data. System (|I]) is referred to as the 
"forward" system and the same system changing t to T — t 
the "backward" system. The indice k introduced below refers 
to the fc*'' back and forth iteration of the algorithm. The letter 
'f stands for "forward" and "b" for "backward", p is the 
designed observer 
For the forward system consider: 



Piit) = ~i[B4t)a, + By{t)ay,pi{t)] 

+ r(a,p{(i)a, -p{(t)) 

- (r + ^)a,iyiit) - y{t)) 

Viit) = Tr(a.p{(f)) 



(2) 



For the backward system consider: 



^A[t) = i[B4T-t)o^ + By{T-t)oy,pl{t)] 



dt 



- r(o-^pt(i)o-z - Pfc(t)) 

+ (r-7V.(yfcW-y(r-i)) 



(3) 



Vkit) 



Tr (T;,/3fe(t) 



Noting that piiO) = p{(T) and p{(0) = pti(r) and 
7>0. 

We initiaHze the observer p^ (0) = /5(0) to be Hermitian 
and of trace 1. We typically take /5(0) = ^I where / is the 
identity matrix. This way we make no a priori assumption 
on the initial state. 

Remark 1: Notice that the observer introduced above is 
trace preserving and stays Hermitian for all time. However 
it does not preserve positivity. 

We define for all t e K+ (fc is defined as: k = E{^) 
where E represents the integer part) 



pit) = piit ~ 2kT) - p{t - 2kT) 

if t e [2kT, {2k + 1)T[ 
pi{t - {2k + l)r) - p{2{k + i)r - t) 

if f e [(2fc + i)T,2(fc + i)r[ 

fixW = B^{t~2kT) 

if t e [2kT, {2k + 1)T[ 
B^{2{k + l)T - t) 

if f G [(2fc + l)r,2(fc + l)r[ 
By{t) = By{t-2kT) 

if t e [2kT, {2k + 1)T[ 

By{2{k + 1)T ~ t) 

if f G [(2fc + l)r,2(fc + l)r[ 
Z{t) = Tr (a.p(i)) 

Let 

V : A (hermitian) -^ Tr [A^) 

V is definite positive. For all A; e N we have: 

^V{p{t)) = -Arv{p{t))--2f{z{t)f 

if t G [2kT, {2k + 1)T[ 
= AVV{p{t))-2^{Z{t)f 

if iG [(2fc + l)T,2(fc + l)T[ 



(5) 



(6) 



(7) 



We note, for all fc e N, 14 = V{p{2kT)). We define the 
function g such that for all t e R+ 



git) 



g4r(t-2feT) .J ^ ^ j2fcr, {2k + 1)T[ 



-4r(t-2(fe+l)T) 



if t G [(2fc+l)r, 2(fc + l)r[ 



Vt. 



fc+l 



/.{2fe+l)T 

Vfc = -27(/ g{t)z\t)dt 

J2kT 
^2(fe+l)T 



+ 



g{t)Z\t)dt) 



(8) 
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Fig. 1. Notice that the upper envelope (i.e Vk = Tr ((p(2fcr))^)) goes 
to zero when k goes to infinity. 
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Fig. 2. Measurement and estimated measurement versus time. 10% 
gaussian noise was added to the data. The estimated measurement is 
obtained by simulating M) with p(0) = pj. (0) and k = 10. 



(Vfe)fe is a decreasing sequence which is studied in more 
detail in Section BM. Before looking into the convergence 
proof, we present some simulations which show the robust- 
ness of the convergence of jO^(O) towards p{0) when k goes 
to infinity. 

III. Simulations 
For the simulations of figures |l| |[ ^ and ^ we take: 

r = 0.25 kHz 7 = 0.25 kHz 

Vi ^B.^{tf + By{ty = 10 kHz T = 1 ms 

We take B^,{t) = Bq cos(6i(t)) and By = Bq siii(6i(t)). So = 
10 kHz. For 9{t), at 10 equally spaced times between and 
T, we take a random value between and 2ti. Using a cubic 
spUne interpolation, we build 6{t) over [0,r]. 

10 iterations are simulated: k = 0,..,10. 10% gaussian 
noise was added to the measurement, B^ and By. 

We initialize the estimator in the completely mixed state 
/5q(0) = p(0) = ^I. We randomly initialize p on the 
Bloch sphere by taking a random Hermitian positive matrix 
satisfying: Tr (p(0)2) = Tr (p(0)) = 1. 
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Fig. 3. The density matrix of the system and its estimator at time t = 0. 
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Fig. 4. The magnetic fields B^it) and By{t). 
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N 



Vo = -27 



(2N+2)T 



g{t)Z^{t)dt 



Vi e M+ g{t) > 1, hence J^ ^ ' Z^t)dt < ^^^. 
Since Vt G R+ ^^(i) > 0, /q°° Z'^{t)dt exists and is finite. 

From (|) and (0) we have Vu G [0, 2T] and Vfc G N: 

V{p{2kT + u)) < V{p{2kT)) (9) 

hence, for all t G M+ we have V{p{t)) < Vq. p is therefore 
bounded and belongs to the ball centered around and of 
radius Vq. 

We are now going to prove that {Vk)k converges to zero 
when k goes to infinity, and from (Q) we will conclude that 
V{p{t)) converges to zero when t goes to infinity. In order 
to prove the convergence of {Vk)k, we are going to prove 
that Z{2kT), f^Z{2kT+), j^Z{2kT+) all converge to zero 
when k goes to infinity. 

We consider B^^^By C'^{[Q,T],M) functions. 



IS 



C^ 



over 



S 



< TD _^ t? _^ R ^ R _^ R 

a^'^y dt-^^^ dt^y^ df^^^^ dt^^y 



[Jkef^]kT,ik + 1)T[. 
and p are bounded 



over S, therefore -^p, -^p, -^p are bounded over S and 
therefore 4:Z,4tiZ,4t^Z are bounded over S. Since Z 

at ' dt^ ' dt-^ 

is continuous over R+ and -^Z is bounded over 5, Z 

dt ' 

is uniformly continuous over IR+, so Z^ is uniformly 
continuous over M+. What's more J„ Z^{t)dt exists and is 
finite. We can conclude by applying Barbalat's lemma [9] 
that limf '^2'' 



Z^{t) =0 and hence 



IV. Convergence proof 

Theorem 1: For any B^.By G C'^{[Q,T],M.) such that 
BM^ByiQ) - By{Q)f^B^{Q) ^ and Vp(0) which is 
Hermitian and of trace 1, we have 

lim Tr(p2(i)) =0 
Remark 2: The convergence outlined in theorem |1| is in 
two steps: 

1) limt^oo Z{t) = 

2) limt^„oTr(p2(t)) =0 

The observer @^ is designed such that we always have 
convergence of the estimated measurement to the measure- 
ment: limt^txj Z{t) — 0. Since the system is observable: 
(cz, CTa;, ay) and its commutators span the space of all trace- 
less 2x2 Hermitian matrices, we can find fields {Bx,By) 
such that limt^oo Z{t) = implies limt_j.oo Tr (p^(<)) = 0. 
Proof: 
For any piecewise continuous function / we define: 

f{2kT+) = \imt^2kT,t>2kT fit). 

From (ph, we know that {Vk)k is a decreasing sequence. 
Besides, for all fc G N, 14- > 0, hence {Vk)k converges to 
a limit that we note by 14c. Summing (||) between and 

iVGN*: 



And in particular 



lim Z{t) = 



lim Z{2kT) = 



(10) 



(11) 



Since the derivatives of Z are not continuous over M+ 
but only over S, we cannot directly apply Barbalat's lemma 
to iZ{t) and £z{t). 

Suppose that jiZ(t) does not converge to zero when t 
goes to infinity. 

There exists e > and a sequence (t„)„ such that 
lim„^oo i« = oo and -^^Zit^) > e (or j-^Z{tn) < -e which 
can be treated in exactly the same way). 
Since -^Z{t) is bounded over S we have: 3?^ G]0,r/2[ 
such that for all n G N and t G [-C™>C"''] lli^i^^ + 
t) - f,Z{t^)\ < |. Where t™" = min(t„ - E{tjT)T,r,) 
and t™"^ = min((^(i„/r) + 1)T- t„, rj). E represents the 
integer part. 

Hence, for all t G [-C™,^''] we have: j-^Z{tn + t) = 
^Z(i„) - ( Az(i„) - Az(i„ + 1)) > ^Z(i„) - I |Z(t„ + 
t) - ^Z{tn)\ >£-§ = §. Also, notice that T > 



j.rmn i -t-max \. ^ 
^n "1" ^n — '/■ 



-Z(t)di| 



> 



rye 



>0 

This is in contradiction with (|10|), we therefore conclude 
that 



lim Z(2fcr) = 
lim X(2kT) = 

lim r(2fcr) = 



This is equivalent to limt 
us to conclude that: 



V{p{2kT)) := 0. (|) enables 
lim Tr(p2(i)) =0 



lim -Zit) = 
t^x dt 



and in particular: 



lim —Z{2kr 

k—^oo dt 



= 



(12) 



(13) 



Suppose that ^Z(2fcT+) does not converge to zero when 
k goes to infinity, there exists e > and a sequence (fc„)„ 
such that lim„^oo ^n = oo and ^Z(2A:„T+) > e, since 
■^Z is bounded over 5, there exists < rj < T such that for 
allng NandO << < 77 |^Z(2/s„r+t)-^Z(2/s„r+)| < 
|. Hence: 



|^Z(2A:„T + 7?)-^Z(2fc„r+)| = 



> 



2fc„T+?7 
2/c„T 

?7e 



di2 



Z{t)dt\ 



>0 
Which contradicts (O). Hence: 



lim — jZ(2A:T+) 
fe-i-oc at 



(14) 



We note 



X(f) = Ti{a,pit)) 
Y{t) = TrKp(f)) 

We recall that from (|llj)(ll3D(|l^' 



lim Z{2kT) = 



lim —Z{2kT+) 

A:— >oo at 
rf2 

lim -—Z{2kT+) 



(15) 



Using (|l]), we find that (|l5|) implies: 

lim Z{2kT) = 



lim B^{2kT)Y{2kT) ■ 
lim —Ba:i2kT)Y{2kT) 

k—^oo dt 



■ By{2kT)X(2kT) 
d 



(16) 



dt 



By{2kT)X(2kT) =0 



Notice that B^{2kT) = B^(0) and By{2kT) = By{Q), the 
same holds for their derivatives. We take B^^By such that 
B,(0)^Sj,(0) - Bj,(0)^B,(0) ^ 0. (Ill) implies: 



V. Extension to the multilevel case 

We now consider a system of total angular momentum F. 
The dimension of the system is d = 2F + 1, and the density 
matrix p{Q) belongs to the set of positive d x d Hermitian 
matrices of trace 1. There are therefore d^ — 1 parameters to 
identify. In order to extend the proof of the two level case 
(d = 2) to the multilevel case, we would need to prove that 
Z{t), 2f.Z{t), .., "! ai_2 Z{t) all converge to zero when t goes 
to infinity. If the system is observable, we would be able to 
conclude that we can find a control such that pit) converges 
to zero. Two complications arise from considering a system 
of higher dimension: 

First, the need to extract information from further derivatives 
of the measurement record systematically reduces the robust- 
ness of the state estimation. One direction of improvement 
would be to design a nonlinear observer which preserves 
positivity. Such an observer is given in [10] for a non 
dissipative system (F = 0). The difficulty is to build an ob- 
server which stays positive even in the backwards dynamics 
which is instable due to the dissipation term in F. Such an 
observer would reduce the size of the admissible /9(0)s given 
a noisy measurement record, the robustness will therefore 
be increased. However, no dramatic improvement should be 
expected since the information on some matrix elements 
of p{Q) are hidden in high derivatives of the measurement 
record. 

Second, although the observability criteria insures the ex- 
istence of a control such that limf^oo Z{t) = implies 
limf^oo p{t) = 0, we don't have any well known method to 
find such a control. The higher the dimension, the harder it is 
to find a control which makes the data y{t) informationally 
complete about the initial state p{Q). 

We now give some simulations which show that our BFN 
protocol still works well for a system of total angular 
momentum F = 1 (d = 'S). Consider the system: 



-p{t) = -i[H{t),p{t)] + rv[o]p{t) 



(17) 



yit) 



TriOpit)) 



Where H is the systems Hamiltonian and I?[0] the Lind- 
blad superoperator. We have: H{t) = gFPB{Bx{t)Fx + 
By(t)Fy) + prFx"^. gp, PB, F and (3 are positive constants, 
Bx , By are the controls and Fx,Fy, Fz are the angular 
momentum operators. O is the observable, and we take 
O = VtFz. V[0]pit) = Op(t)Ot_ i(ootp(i)+p(i)oto). 
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Fig. 5. The upper envelope is Vk = Tr ((p(2fcT))2). Notice tliat it 
decreases and seems to converge to zero. 



The superscript ^ stands for conjugate transpose. The term 
fiTFr^'^ is necessary to insure the observabiHty of the system 
[11]. We now consider the following observers: 






m 



= -z[H{i)rpi{i)]+TV[0]pi{t) 



TrfOp{(i) 



y{t)) 



(18) 




Fig. 6. Measurement and estimated measurement versus time. 10% 
gaussian noise was added to the data. The estimated measurement is 
obtained by simulating (|l7|) with p(0) = p^ (0) and k = 50. 



'(0)1 



\P5 



■M 





dt 



pi{t) = i[H{T-t),pl{t)]-TV[0]pl{t) 

- lO{yl{t)-y{T-t)) 
yiit) = TT{Opi{t)) 



(19) 



with the conditions: p^(0) = pliT) and p{(0) = 
/5|_i(T). We initialize the observer at /5q(0) = ^Id where 
Id is the 3x3 identity matrix. 

For the numerical simulations in figures g, |^ and ^ we 
take: 

gp = 1 fJ-sBo = 30 r = 1 

7 = 1 /3 = 10 T =1 

We take B^{t) = Bocos{d{t)) and By{t) = Bosm{9{t)). 
9{t) is found using a numerical search routine aiming to max- 
imize a certain criteria (entropy), as explained in [12]. 10% 
noise is added to the controls B^, i?„ and 10% noise is added 

/I 1^ 
to the measurement y{t). We take p{Q) = | 

VI 1, 

Notice that the estimated measurement is almost identical 
to the measurement y{t) (figure ^. Also, the sequence 
{Vk)k decreases and seems to converge to zero (figure 
Bh. This enables us to reconstruct the initial state with a 
96% fidelity where the fidelity T is computed as follows 



J" = Tr 




pLmp{Q)^piM 



[13]. More iterations 



are needed than in the 2 level case (50 as opposed to 10) to 



Fig. 7. The density matrix of the system and its estimator at time t = 0. 
We plot the modulus of each matrix element. 



achieve a similar fidelity. Each back and forth iteration takes 
about 0.1 seconds so the presented simulation takes about 
5 seconds to run. As mentioned above, the fidelity of the 
reconstruction can be increased if a better control is found. 

VI. Conclusion 

In this paper we propose a BFN scheme to estimate the 
initial state of a quantum system when a continuous mea- 
surement of a single observable is given over a time interval 
[0,T]. A convergence proof and simulations are given for 
the two level case, and the considered experimental settings 
were similar to those in [5]. We discuss the extension of 
this algorithm to the multilevel case outlining the limitations 
and possible improvements of this protocol, and we present 
simulations in the case of a spin 1 system. A quantitative 
comparison of this method to the ones considered in [4] 
and [6] in terms of estimation time and robustness will 
be necessary to put forward the advantages of this state 
reconstruction protocol. 
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